Molecular dynamics simulations of spin and pure liquids with preserving all the 

conservation laws 
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A new methodology is developed to integrate numerically the equations of motion for classical 
many-body systems in molecular dynamics simulations. Its distinguishable feature is the possibility 
to preserve, independently on the size of the time step, all the conservation laws inherent in the 
description without breaking the time reversibility. As a result, an implicit second-order algorithm 
is derived and applied to pure as well as spin liquids for which the dynamics is characterized by the 
conservation of total energy, linear and angular momenta as well as magnetization and individual spin 
lengths. It is demonstrated on the basis of Lennard- Jones and Heisenberg fluid models that when 
such quantities as energy and magnetization must be conserved perfectly, the new algorithm turns 
out to be more efficient than popular decomposition integrators and standard predictor-corrector 
schemes. 
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I. INTRODUCTION 



During the last years, considerable attention has been 
focussed on computer experiment studies of relaxation 
properties and critical phenomena in classical spin sys- 
tems . These studies dealt mainly with lattice mod- 
els such as the Ising, the XY, and the Heisenberg. Of 
notable current interest is the investigation of continuum 
spin liquid models ]To|-p^ in which additional dynamical 
effects are possible because of the coupling between spin 
and liquid subsystems. 

Quite recently |16 1^, a set of symplectic algorithms of 
different orders in the time step has been constructed for 
numerical integration of motion at the presence of both 
translational and spin degrees of freedom. As a conse- 
quence, the molecular dynamics simulations (MD) of a 
Heisenberg spin fluid have been carried out for the first 
time. The symplectic integrators were derived by devel- 
oping the Suzuki- Trotter technique for decomposi- 
tions of exponential operators. Their main advantages 
over standard predictor-corrector schemes are explicit- 
time reversibility and exact conservation of spin 



ness, 

lengths. It was also shown that the decomposition al- 
gorithms permit significantly larger time steps and lead 
to a substantial speedup of the calculations. In a par- 
ticular case when the spin degrees of freedom are frozen, 
these algorithms can be reduced to the well-known ve- 
locity Verlet integrator p9[ |, widely used for simulating 
of pure liquid dynamics. 

However, the decomposition approach (like predictor- 
corrector and other existing traditional numerical 
schemes such as Runge-Kutta, etc.) do not pre- 

serve the total energy and magnetization of the system. 
In most MD applications the accuracy achieved for the 



energy-magnetization conservation by the decomposition 
algorithms is high enough in order to obtain reliable re- 
sults. Moreover, this accuracy can be improved using 
higher-order versions of the decomposition approach 
or decreasing the step size. But if the integrals of mo- 
tion have to be conserved perfectly, the non-conservative 
algorithms may not be an optimal choice for the solu- 
tion of the problem. The reason is that then the time 
step should be divided into a lot of subintervals, reduc- 
ing considerably the efficiency of the computations. 

The exact conservation of integrals of motion is espe- 
cially important for simulations of spin liquids near phase 
transitions, when the phase diagrams, dynamical scaling, 
long-time correlated behavior or derivatives of the ther- 
modynamic functions are investigated. This is so because 
in these cases, the presence of artificial fluctuations in en- 
ergy and magnetization may have a signiflcant influence 
on the results. Therefore, it is desirable to look for an 
algorithm which conserves the fundamental physical in- 
variants exactly or at least within machine accuracy. 

In the present paper, a novel approach to numerical 
integration of the equations of motion for spin and pure 
liquids is introduced. The main feature of this approach 
is its intrinsic preservation of all the conservation laws 
inherent in the system without violating the time re- 
versibility property. The paper is organized as follows. 
The basic equations and their integrals of motion are 
described in Section 2. Section 3 is devoted to a conse- 
quent derivation of the desired second-order algorithm. 
Its possibility to conserve exactly the integrals of motion 
is demonstrated there also. In Sections 4 and 5, the algo- 
rithm is tested in actual MD simulations on Heisenberg 
and Lennard- Jones fluid models, respectively, and com- 
pared with previous numerical schemes. The discussion 
and concluding remarks are given in Section 6. 
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II. BASIC EQUATIONS OF MOTION AND 
CONSERVATION LAWS 



Let us consider a classical A^-body system described 
by the Hamiltonian 



H 



N 



i=l 



1 ^ 



where and Vi are the translational position and ve- 
locity, respectively, of particle i with mass rrii carry- 
ing spin Si. The fluid part of the potential is denoted 
by ip{rij), whereas J{rij) is the exchange integral cor- 
responding to a pair of spins with the interparticle sep- 
aration rij = \Ti — Yj\. Note that within the classical 
approach, each spin is treated as a continuous three- 
component vector with fixed length (putting for conve- 
nience I Si I — 1, so that J will be measured in energy 
units). Although the results which will be obtained below 
can easily be adapted to a larger class of Hamiltonians 
(to multi-component systems, for instance), we restrict 
ourselves for the sake of simplicity to the basic model 
(1) which represents a typical isotropic Heisenberg spin 
fluid dmil]. For J = 0, Eq. (1) reduces to a pure hquid 
model. 

In MD simulations it is necessary to solve numerically 
the equations of motion, dp/dt ~ [p,H], where [, ] de- 
notes the Poisson bracket and p = {ri,Vi,Si} is the full 
set of microscopic phase variables. For the system under 
consideration, the dynamical equations can be written 
more explicitly p^ , p7[ , 

dr,- 



dt 



dvj _ Ji_ ^ _ J_ / dip, 
dt rrij rrij ^-^ ^ \ dri 



^"^'^s,.s,)i^, (2) 



-XK. = -X ^ J,, 



\ dri J dn^ 



dt h 



Here, fi 



^JiJ^^) and g, 



E 



gy are the force 



and internal magnetic field, respectively, acting on parti- 



cle i due to the interactions U 



-(<PL-'/I,s»-Sj)ry/r,. 



and gij = Jij Sj with all the rest of bodies, where 
ifij = f{rij) and Jy = J{rij). Note that the quantum 
Poisson bracket was applied [p|Jl^ to derive the equations 
for spin subdynamics. If an initial state p(0) is specified, 
the time evolution p{t) can be uniquely obtained by in- 
tegrating Eq. (2). 

Taking into account the symmetry ipij = ipji and 
Jij — Jji of interaction potentials, it follows from Eq. (2) 
that the total energy E = H, the total magnetization 
M = J2i the total linear momentum P = J2i n^i'^i as 
well as angular momentum L = mir^ Xv^ are integrals 
of motion, i.e., dE/dt = dM/dt = dP/dt = dL/dt = 0. 
The structure of spin equations of motion (the last line 



of Eq. (2)) imposes in addition the conservation of indi- 
vidual spin lengths, |si| = Si = const. Indeed, dsi/dt = 
d(si-Sj)i/2/dt = {s.rdsi/dt)/si = s^-[siXgi]/h.Si = 0, be- 
cause the equality a'[axb] — is valid for arbitrary 
vectors a and b. Besides, the exact solutions are time 
reversible, since the equations of motion arc invariant 
with respect to the time inversion transformation t — > — t, 

{V^,sj {-Vi, -sJ. 

No existing numerical scheme can obey perfectly all 
the just mentioned properties. The exact conservation 
during the integration can be achieved only for some of 
the integrals of motion, such as linear momentum, for 
example. Usually, it is required that the deviations of 
conservative quantities from their original values to be 
within an acceptable level of precision. This results, how- 
ever, in limitations on the size of the time steps which 
actually can be used for MD simulating. 



III. THE METHOD OF INTEGRATION 

We will show in this section that it is possible to 
generate time-reversible microscopic trajectories along of 
which all the integrals of motion are preserved at arbi- 
trary finite time steps. Our derivation of the desired algo- 
rithm is started by considering a mid-point scheme of the 
second order. According to this scheme, the dynamical 
variables can be propagated as 

p{t + r) = p{t) + T[dp/di]i+,/2 + O(t') 

with T being the step size and 0{t^) denoting the trun- 
cation terms. In view of Eq. (2), the explicit expressions 
for such a propagation read 



r,(t + r) =r,(i) + rv,(t+§), 
v,(t + r) = v,(t) + Tf,{t + §)/m,, 
Si{t + r) = Si{t) + T[siXgi]t+^/U, 



(3) 



where the mid-step values of , and s.; X gi should be 
specified additionally. 

The only way to construct mid-point translational ve- 
locities maintaining the time reversibility property is 

TV,{t + i) = ^ [v.(i) +Mt + r)] + 0{t^). (4) 

The terms 0{t^) of third- and higher-orders can be ig- 
nored because the corresponding terms of the same or- 
ders have been truncated already within the mid-point 
approach. Eq. (4) represents, in fact, an implicit inter- 
polation formula in which past (at time t), and future (at 
t + r) values of dynamical quantities enter symmetrically, 
assuring automatically the reversibility of the solutions. 

In the case of translational forces, there are several 
possibilities to build the mid-point values. The reason 
is that the interparticle function iij = f (r^j, s^-Sj) de- 
pends explicitly on relative position r^j and orientation 
Si'Sj which in turn vary with time. Thus, we can apply 



2 



the mid-point interpolation either to the function fy as 
a whole, or directly to the dynamical variables rj^ and 
Si'Sj. As a result, two different approaches to the force 
evaluation may be introduced, namely, 

fij(*+i) = \ f [Si-Sj]t) +f{rij{t + T), [Si-Sj]t+r] 

and 

%(i+i)=f(rii(i+i),[si-s,]t+x). (5) 

The last approach requires the knowledge of mid-point 
values for ani 
tive positions is 



values for and Sj-Sj. The obvious choice for the rela- 



(6) 



The interpolation of the scalar product [sj'Sj]t_|_^/2 will 
be described latter. None of the above approaches for 
evaluating fij{t + r/2) can lead to a scheme with exact 
preservation of the total energy. The energy will only be 
conserved approximately with the precision within which 
the microscopic solutions are calculated, i.e., E{t + r) = 
E{t) + 0{T^). 

We will show now that the second approach (Eq. (5)) 
can be modified in such a way as to compensate the loss 
of precision in the total energy. The idea lies in the fol- 
lowing. Since, according to Eq. (1), the energy difference 
E{t + T) — E{t) is a function of the quantities (p{rij (t + r)) 
and ^{Tij{t)) as well as J{rij{t + r)) and J{rij{t)), it is 
natural to try to evaluate numerically the partial deriva- 
tives ip'{rij) = d(p/drij and J'{rij) = dJ/drij (which 
appear in Eq. (5) at t + T/Q.) in terms of the same quanti- 
ties, rather than to calculate the derivatives analytically. 
This is possible because for any fimction £,{rij) depend- 
ing only on the interparticle distance Vij we can write the 
following two expressions 



d^(rij) 5^ drj 



di 



and 



^{nj{t + T))-£,{n,{t)) 



combining of which gives 



^{rij{t + r))-£,{n,{t)) 



+ 0{t% (7) 



where the mid-point values of relative velocity 
Vj — Vj are calculated according to Eq. (4) as 



Then choosing ^ = (p,J, one finds the expression 



rii(i+i)-Vi,(i+i) 



-'fiinAt)) - [J{rij{t + r)) - J{rij{t))] [srs,]t+x j (8) 

for mid-step translational forces, where the 0{t^) terms 
have been neglected. 

Performing scalar multiplication of Eq. (8) with the 
vector Vj(f -|- t/2), then taking the sum over all the par- 
ticles [i = 1,2, ... ,N), and using the fact that the double 
sum obtained in the right-hand side is invariant with re- 
spect to the replacements z j, it can be shown that 

N 1 ^ / 

f.(t + i).v,(t + I) = ^ ^{n^{t + r)) 



i=l ijij 

-finjit)) - [J{rij{t + T)) - J{rij{t))][si-Sj]t+i 



Assuming for the moment that spin degrees of freedom 
are frozen (i.e., that [Sj-Sj] do not depend on time), the 
last relation can be presented in the form 

N 

T ^ fi(i + T/2).Vi(i + t/2) = U{t) - U{t + t), 

i=l 

where U denotes the potential energy of the system. On 
the other hand, multiplying the second line of Eq. (3) by 
Vi{t + t/2) and summating over the particles yields 

E ^ + ^) - it)) • (Vi {t + T)+ Vi (i) ) 

i 

N 

^ K{t + r) - K{t) =rJ2 ^ii* + i)-Vi(i + i), 



i=l 

where K denotes the kinetic energy. We sec, therefore, 
that during the time propagation given by Eqs. (3) and 
(8), the total energy E = K + U is conserved exactly for 
any t, i.e., E{t + t) = E{t). 

Note that Eqs. (7) and (8) are well defined when 
the scalar product Yij{t -\- r/2)-Vy(t -|- r/2) tends to 
zero. This is so because according to the first line of 
Eq. (3), the mid-step relative velocity is connected with 
the change in position by the constraint Wij{t + t/2) = 
{vij{t+T)—Vij{t))/T. So that the scalar product is merely 
equal to {Yij{t -\-t)-\- Vij{t))-{Yij{t -\-t) - Vij{t))/{2T) = 
{nj{t + t) + r,j{t)){r,j{t + r) - (t))/(2T). As a resuh, 
the right-hand side of Eq. (7) can be rewritten in the 
following mathematically equivalent form 



i{nj{t + T))-^{nj{t)) 



Tij (t + T)- nj (t) Uj (t) + Tij [t + t) 



+ 0{t^) (9) 



which reduces to ^'(r^ (< + T/2))/rij{t + r/2) + e'^0{T'^) 
when the value \rij{t + t) — rij{t)\ < e is small enough, 
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where rij{t + T/2) = {rij{t)+rij{t + T))/2 and denotes 
a machine zero. 

The exact energy conservation can also be achieved in 
the presence of spin subdynamics. In order to show this, 
unfreeze now the spin variables, and consider first the 
question of how to interpolate the vector product SjXgi 
arising in the third line of Eq. (3). Again, since this prod- 
uct depends on time implicitly via dynamical variables 
and gi, we will have here a lot of possibilities. The first 
of them 

[SiXgi]t+x = i Si{t)Xgi{t)+Si{t + T)Xgi{t + T) 

is not suitable because it does not lead to the con- 
servation of individual spin lengths, i.e., Si(t + t) = 
Si{t) + 0{T^). At the same time, the second interpolation 

[SiXgi]t+x = i[s,(i)+Si(t + T)lxgi(t+f) (10) 



does conserve spin lengths exactly, Si{t + r) = Si{t), 
for arbitrary choice of gi{t + r/2). Indeed, substitut- 
ing Eq. (10) into the propagation equation Si{t + r) = 
Si{t) + r[siXgi]t+x/2/?i and solving analytically the ob- 
tained expression with respect to Si{t + r), one obtains 



Si(i + T) 



1 + 



1^ 

it+i 



^(2 . ([g,],+5.s,(t)) - )%(t)))]. (11) 

As can be verified readily, Eq. (11) represents an unitary 
transformation, Si{t + t) = &i{t,T)si{t), where &i{t,T) 
is a rotation matrix which, of course, does not change the 
norm of vectors. 

Three different time-reversible interpolations can be 
introduced for the factor [gi]t+r/2 = Ei{t + 't/'^) = 
T,jUjti)\.Sij]t+r/2- They are: 



J{n,{t))sj{t) + J{nj{t + T))sj{t + t) 



^ ( nj jt) + nj (t + r ) Sj jt) + Sj {t + t) 
[Sijlt+i — J y ^ j ^ , 



and 



[Eio]t+i 



J {nj {t))+J {nj {t + t)) Sj jt) + s,- {t + r) 



Eqs. (3)) over all the particles and taking into account 
Eq. (10) gives M{t + t) = M{t) + AM, where 

= ^ E Jij{t+i)h{t)+Si{t+r)] X [Sj{t)+Sj{t+T)] 



and Jij (t + t/2) may be equal either to J{{rij (t) + Vij [t + 
t))/2) or (Jinjit)) + Jin-iit + t)))/2. The term AM is 
canceled because of the invariance of the double sum with 
respect to the transformation i j, and of the obvious 
equality axb-|-bxa = which fulfills for any vectors a 
and b. Thus, M{t + t) ^ M{t) in both the cases. How- 
ever, in the first of them when Jij{t + r/2) = J{{rij{t) + 
rij{t + T))/2), the energy difference E{t + T) — E{t), being 
a function of two quantities J{rij{t + r)) and J{rij{t)), 
cannot be reduced to zero exactly using only one mid- 
point value J{{rij{t) + rij{t + t))/2). 

At the same time, within the last third interpolation 
given by Eq. (12) we are able to perform such a reduction. 
To demonstrate this, let us consider finally the interpo- 
lation of the spin scalar product [si'Sj]t+r/2 appearing in 
Eq. (8). Similarly to Eq. (6), one defines this interpola- 
tion in the form 



Si{t)-Sj{t)+S,{t + T)-Sj{t + T) . (13) 



Then, using Eqs. (3), (8), (10) and (12), it can be shown 
that the following equality holds 



N 



N 



where U^^^ denotes the spin part of the potential energy. 
So that, as in the case of frozen spin subdynamics, the 
sum Tj2iLi + T/2)--Vi{t + t/2) is reduced to the po- 
tential energy difference U{t) — U{t + r). But as was 
shown earlier using Eq. (3), this sum can be expressed 
also as the difference K{t+T)—K{t) in the kinetic energy. 
This indicates again that the total energy E = K + U is 
conserved exactly, i.e., E{t + t) = E{t), despite the mi- 
croscopic solutions ri{t + t), Vi{t + r) as well as Si{t + r) 
are obtained with a limited 0{t^) accuracy. It is worth 
mentioning that the interpolation 



(12) 



Si (t) +Si{t + T) Sj (t) + Sj {t + t) 



The first approximation cannot be chosen for our purpose 
because it destroys the total magnetization of the system, 
i.e., M(t+T) = M(t)+0(T^). The last two interpolations 
do reproduce the magnetization vector perfectly. Indeed, 
summing up the individual spin propagation (third of 



instead of Eq. (13) is possible, in principle, too but it will 
not lead to the energy conservation. 

The approach considered conserves also the total lin- 
ear and angular momenta, i.e., P(i + t) = P(t) and 
'L{t + t) = L(f). The first follows directly from the 
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structure of mid-point translational forces (8) for which 
X;,f»(i + T/2) =0, so that 

rriiVt + t) = ^ niiVi (t) . 



Further, taking into account Eqs. (3) and (4), the posi- 
tion propagation can be cast as 

r,{t + t)= Y,{t) + v,(t)r + f,(t + f )r72m,. 

Then the sum + T)xvi{t + r) reduces to 

m^r,{t)xvi{t)+T r^{t+T/2)x{^{t+T/2). The last 
term is canceled since, according to Eq. (8), the interpar- 
ticle forces are parallel to mid-step vectors rij{t + t/2), 
and, thus, the second property 

m,Y, (t + r) X V, (i + r) = ^ m,r, (i) X [t) 



is also satisfied. 

Thus, the desired algorithm of the second order has 
been constructed. In view of Eqs. (3), (4), (6)-(10), (12) 
and (13), the algorithm can be presented in the following 
compact form 



v,{t + t)^ v,(t) + ^ [v,(<) + v,(< + t) 
v,(t + r) =v,(t) 



T_ r,j-(t) +r,j-(t + r) ^ 
ip{r,j{t + r)) - (^(r„ (t)) J(ry (t + r)) - J(r„ (i)) 



ry(t + r) -ry(t) 
Sj(i)-Sj(t) + Sj(t + T)-Sj(t + r) 



At) 



(14) 



s.(t + r)^s.(t) + I!ii^)±|^X 



E 



J{nj{t)) + J{nj{t + t)) Sj{t) + Sj{t + t) 



Equation (14) constitutes, in fact, a coupled system of 
three nonlinear vector equations for each particle with re- 
spect to the same number of unknowns ri(t-|-T), 'Vi{t + T) 
and Si {t + t). The system can be solved in a quite effi- 
cient way by iteration, letting initially v^°^ (* + t) = (t) 

and sf'\t + T) = Si{t) in the right-hand sides of Eq. (14). 
Then the current values for ri{t+T), Vi(/;+r) and Si{t+T) 
obtained in the left-hand sides of Eq. (14) are treated as 
initial guesses for the next iteration. Already two it- 
erations are sufficient to reach the 0{t^) accuracy for 
the microscopic solutions and energy conservation, i.e., 
Eit + r) — E{t) — 0{t'^). The goal of carrying out further 
several updates of Eq. (14) is to reduce the uncertainty 
e = E{t + t) — E{t) in energy deviation to a negligibly 
small value (by adjusting the number Z > 2 of iterations 



for a given r). The rapid convergence e —>■ +0 is guaran- 
teed by the relative smallness of the step size r and an 
exponential decaying of e with increasing / (see the next 
section). 

It is interesting to remark that the total linear and 
angular momenta are conserved exactly within each it- 
eration of Eq. (14). The reason is that the interparti- 
cle forces are evaluated exploiting Newton's third law 
and the velocities Vi{t + r) are updated before all (i — 
1,2,..., TV) the advanced positions ri{t + t) were cal- 
culated. For similar reasons, the magnetization con- 
servation is also fulfilled for each iteration, when the 
spins are updated according to the third line of Eq. (14). 
In this case, however, the individual spin lengths will 
only be preserved like energy in an iterative sense, i.e., 
Si{t + t) — Si{t) = 0{e). Nevertheless, the spin lengths 
can be maintained exactly within each iteration by re- 
placing this third line by its mathematically equivalent 
counterpart (11) (were [gi]t+T/2 being evaluated with the 
help of Eq. (12)). 



IV. NUMERICAL TESTS. COMPARISON WITH 
OTHER METHODS 

In our MD simulations of the Heisenberg spin fluid 
(Eq. (1)), the Yukawa function III 



y \ r) ~ w— exp 

r \ (J 



(15) 



was used to describe the spin-spin interactions. Thcji 
uid subsystem was modelled by a soft-core potential 



(/5(r) — Au 



12 



(16) 



which accepts nonzero values at r < 2^/^cr and Lp{r) = 
at r > Here cr is the diameter of particles, and u as 

well as w denote the intensities of core-core and spin-spin 
interactions, respectively. The simulations were carried 
out in a microcanonical ensemble for N — 1000 identical 
{mi = m, Si = 1) particles in a cubic box of volume V = 
employing periodic boundary conditions. The Yukawa 
function was truncated at = 2.5f7 < L/2 and shifted 
to be zero at the truncation point to avoid the force sin- 
gularities, i.e., J(r) = Y{r) — Y{Rc) at r < i?c and 
J(r) = otherwise. We have chosen the same thermo- 
dynamic point as considered in previous papers |]l6| , |l7| , 
namely, a reduced density of n* = Na'^ /V — 0.6, a re- 
duced temperature of T* = k^T/w = 1.5 < T* (where 
k-Q is the Boltzmann's constant and T* « 2.055 the crit- 
ical temperature of the system ^^ ), a non-zero magne- 
tization per particle of |M|/A^ = 0.6536 as well as the 
same values for the reduced core intensity u/w = 1 and 
the dynamical coupling parameter d = a{m'wy^'^ /h — 2. 

The equations of motion were solved using a well es- 
tablished Adams-Bashforth-Moulton (ABM) predictor- 
corrector integrator of the fourth order pffl, the explicit 
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decomposition schemes [^jl? of the second (ED) and 
forth (ED4) orders, as well as our conservative spin fluid 
dynamics (CSFD) algorithm (Eq. (14)). All the test runs 
were started from an identical well equilibrated conflgu- 
ration. A typical example for the reduced total energy 
E* = E/w and magnetization M = |M| per particle as 
depending on the length of the simulation is shown in 
subsets (a) and (b) of Fig. 1, respectively, at a reduced 
time step of r* = r(w/m(T^)^/^ — 0.01. 

The huge systematic deviations in the total energy ob- 
tained within the ABM approach (see the dashed curve 
in Fig. 1 (a)) points out clearly that it is highly unsta- 
ble and, thus, not suitable for long-duration observations 
over the system at the time step considered. We men- 
tion that in the ABM scheme, the dynamical variables 
are first predicted, 

p{t -f t) = p{t) + [55p(t) - 59p(t - t) 

+ 37p(i - 2t) - 9p(t - 3r)] ^ + 0{t-), 

and further iteratively corrected as 

p{t + t) = p{t) + [9p(t + t) + 19p(i) 



bp{t -t) + p{t - 2r) 



24 



where p — {vi,ii/ni, [siXgi]/fi}. The strong instability 
of the ABM integrator can be explained by the facts that 
it destroys the unit norm of spin lengths (although con- 
serves the magnetization vector) and generates time ir- 
reversible solutions (as is now rigorously proved , the 
numerical stability follows directly from the reversibility 
of an algorithm). For this reason, the ABM as well as 
other existing predictor-corrector schemes can be used 
only at very small time steps (namely, at r* < 0.00125, 
see Refs. ||l^,|l^), were they exhibit similar equivalence 
in the energy conservation as that of the decomposition 
algorithms. However, such small step sizes are inefficient, 
because then too much expensive force and field recalcu- 
lations have to be performed in order to cover the fixed 
observation time. Moreover, the ABM scheme is about 
twice times slower compared to the ED integrator even 
if one iteration only is applied within the corrector pro- 
cedure. 

No drift in the functions E{t) and M{t) was recognized 
within both the decomposition ED and ED4 algorithms 
at time steps up to r* = t*^,^^ = 0.01 and over a num- 
ber of time steps of t/r = 100 000. These algorithms, 
however, do not conserve the total energy and magneti- 
zation exactly. Instead, the last functions fluctuate quite 
visibly especially in the ED case, as can be seen from 
Fig. 1. The ED4 energy fluctuations are approximately 
in factor two smaller than those of the ED algorithm. 
This compensates to some extent the additional proces- 
sor time needed to evaluate high-order expressions. But 
the ED4 algorithm allows to reduce the magnetization 
deviations to a negligible small level which does not ex- 
ceed about {{M{t) - M(0))2)i/2/7v - 3 • 10"^ even at 



the greatest time step rj^j^xj where ( ) denotes the mi- 
crocanonical averaging. It is worth mentioning that the 
ED/ED4 energy-magnetization fluctuations are caused 
by the 0{t^)/0{t^) truncation errors and thus they will 
increase drastically with increasing r. 
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FIG. 1. The total energy E* /N (subset (a)) and magnetiza- 
tion M/N (subset (b)) per particle as functions of the length 
t/r of the simulations performed for a Heisenberg spin fluid 
using the predictor-corrector (dashed curve in (a), marked as 
ABM), decomposition (solid curves, ED/ED4), and our new 
(bold solid horizontal lines, CSFD) algorithms. Note that the 
ABM and ED4 curves are indistinguishable in (b) from the 
CSFD line. 

The situation is completely different in the case of our 
new approach, because the CSFD algorithm preserves the 
integrals of motion for arbitrary time steps. Of course, 
we cannot apply too large step sizes {t* ^ 1) since then 
the microscopic solutions will deviate considerably from 
their exact counterparts and because of too large number 
of iterations needed to achieve the convergence. Choos- 
ing, for instance, t* = r^^x = 0-01 we have determined 
the following levels £ in the averaged total energy fluc- 
tuations £ = {{E*{t) - E*{O jfy/^/N at the end of the 
100 000 time step runs: 9.2 • 10"'*, 2.3 • 10-^ 3.1 • 10-^ 
2.2-10^^, and 2.8 -10^^ corresponding to the numbers I of 
iterations: 2, 3, 4, 6, and 8, respectively. We see, there- 
fore, that the iterations converge rapidly with increasing 
I and the uncertainties can be approximately described 
by the exponential dependence f ~ 3 • 10"'' exp(— 1.2Z) 
at / > 4. Of course, the iterative solutions require addi- 
tional computational efforts, but they are justified when 
a high level of the energy conservation is necessary. In 
order to demonstrate this, we have tried to reduce the 
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energy fluctuations within the ED/ED4 algorithms by 
decreasing the time step. The corresponding result for E 
at the time steps 0.01, 0.005, 0.0025, and 0.00125, i.e., at 
r* = T^^^/l with / = 1, 2, 4, and 8 is presented in Fig. 2 
in comparison with the above CSFD data. 
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FIG. 2. The averaged total energy fluctuations f as a func- 
tion of the number I of iterations, obtained in the Heisenberg 
fluid simulations within the CSFD algorithm at the reduced 
time stop Tniax = 0.01 (bold solid curve). The levels of £ cor- 
responding to the decomposition integrators at the time steps 
r* = T^ax/l are connected by solid (ED) and dashed (ED4) 
curves. 

As one can see, such a reduction of the ED/ED4 en- 
ergy fluctuations is not efflcient, since the deviations £ 
behave like ~ i.e., decrease with increasing I much 
more slower than the exponential dependence obtained 
within the CSFD algorithm. In view of the results of 
Fig. 2 and taking into account that at the same value of 
I one needs approximately the same processor time with 
both the ED and CSFD algorithms (the ED4 integrator 
needs the time larger by factor 5 and is less economi- 
cal) in order to investigate the system over an identical 
time interval, we come to the following conclusion: When 
the total energy must be conserved up to a precision of 
£q ~ lO^'' (the intersection point of the ED and CSFD 
curves, see the horizontal dashed line in Fig. 2) or better, 
the preference should be done to the CSFD algorithm. 
For example, a level of f ~ 10~^ in the conservation is 
achieved at Z ~ 5 within the CSFD algorithm, while up 
at / 50 for the ED scheme (the last value was obtained 
by extrapolating the ~ dependence to larger /). Thus 
the CSFD algorithm appears to be approximately in 10 
times faster than the ED integrator at this level of £. 
For £ > £o,we can restrict ourselves to the usual explicit 
decomposition integrators. 



Note that despite the uncertainties £o 10~^ look 
quite small, they can influence considerably on some 
observable macroscopic quantities. The influence can 
be estimated quantitatively in terms of the ratio F = 
{{{E{t) - E{Q) f)/{{U{t) - C/(0))2))i/2 of total and po- 
tential energy fluctuations. For our system U = {{U* (t) — 
C/*(0))2)i/2/7v ^ 10-2, where U* = U/w, so that 
To = £o/U ~ 1%. Usual investigated quantities, such as 
thermodynamic hmctions, structure factors, etc., will be 
calculated approximately with the same relative precision 
Fo (provided the averaging over the produced trajecto- 
ries is performed during a sufficiently large time interval 
to be entitled to ignore the statistical noise). But when 
long-tail time correlation functions or derivatives of the 
thermodynamic functions are involved into the computa- 
tions, the impact of the artificial energy fluctuations on 
the results will be much greater. For instance, the rela- 
tive uncertainty in the measurements of the specific heat 
(which are based on a microcanonical ensemble fluctua- 
tion formifla) is estimated to be already (ro)^/^ ^ 10%. 
This uncertainty may appear to be too large to determine 
correctly a phase diagram of the system. 

A similar pattern to that shown in Fig. 2 was ob- 
served within the CSFD approach at greater time steps 
r > 0.01. The energy as well as magnetization fluctu- 
ations continued to damp exponentially with increasing 
I, although a greater number of the iterations was nec- 
essary to reach the same level of the conservation. Note 
that a rotational matrix version of the CSFD algorithm 
(when the third line of Eq. (14) is replaced by Eq. (11)), 
in which individual spin lengths are maintained exactly 
within each iteration, leads to a somewhat better energy 
preservation at a given I (but then the total magnetiza- 
tion like energy will be conserved in the iterative sense, 
i.e., at sufficiently large I). The CSFD results presented 
above for £ have been obtained using this version for spin 
subdynamics propagations. 

Further improvements in the efficiency of the CSFD 
algorithm can be reached applying the following com- 
putational trick. It can occur that after a some period 
of time during the integration process the energy differ- 
ence E{t) — E{0) corresponding to the last l-th iteration 
(within a current time step t/r) exceeds the difference 
obtained for the previous (l — l)-th iteration. Such a 
situation is possible because of romid-off errors and an 
accumulation of other numerical uncertainties, especially 
at relatively small values of I, where the lack in the time- 
reversibility can lead to an instability of the solutions 
(note the CSFD algorithm is time-reversible in the iter- 
ative sense, i.e., at large enough values of Z). Then to 
avoid the accumulation, we should take merely the val- 
ues for microscopic phase variables corresponding to this 
previous (I — l)-th iteration. The trick with a flexible 
number of the iterations will guarantee a good stability 
for small I ~ 2-4 as well. 

Another technical detail concerns the way in which the 
expression [^{r{t + T)) — ^{r{t))]/[r{t + T) — r{t)] (appear- 
ing in Eq. (14) for ^ = J) should be treated in the 
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limit r{t + r) — * r(t). As was pointed out earlier, this 
expression must be computed using its limiting represen- 
tation ^'{[r{t) + r{t + T)]/2) + 0{t^ ) when the difference 
\r{t + t) — r{t)\ < e is small enough. Then letting be- 
ing equal a machine zero, the truncated term e^C'(T^) 
can be ignored completely. In our program code we have 
used a double precision throughout with 16 significant 
digits, = 10~^^, so that the value e was set to be equal 
to 10~*. It is interesting to remark that the condition 



, (t)\ < e was never achieved for any pair ij 



of particles during the simulations and, thus, the limiting 
expression was never used. This can be explained by the 
fact that the probability of finding the system in such 
a state is prohibitively small and is proportional to Ce. 
The coefficient C increases with increasing the length of 
the simulations and the number of particles as C ~ tN'^ . 
Thus, the limiting expression is expected to be applied 
for systems with a greater size or when extra long simu- 
lations are performed. 

Finally, some words about the angular momentum con- 
servation. As is well known, the periodic boundary con- 
ditions, which are commonly used in MD simulations 
to reduce the finite-size effects, destroy the angular mo- 
mentum vector. Nevertheless, it has been established 
that this vector is conserved in our simulations in mean, 
namely, (L(i)) « L(0). Note that initial values for the 
total angular as well linear momenta were putted to be 
equal to zero, L(0) = and P(0) = 0, i.e., the system 
was considered at the very beginning as one which does 
not move as a whole translationally and rotationally. 



V. APPLICATIONS TO OTHER SYSTEMS: 
PURE LIQUIDS AND HARMONIC OSCILLATOR 

The algorithm derived in the preceding section can also 
be applied with equal successes to dynamics simulations 
of other liquid models. For instance, letting formally 
J = 0, we come to the usual equations of motion corre- 
sponding to a pure liquid system. These equations can 
be integrated using the first two lines of the same prop- 
agation equation (14), where the terms with J in the 
right-hand side of the second line must be omitted (the 
third line describes spin subdynamics and is not relevant 
in this case). 

Our simulations of pure liquid dynamics were based 
on a system composed of iV = 256 particles inter- 
acting through a cut-off Lennard- Jones (LJ) potential, 
(^(r) = $(r) - $(i?c) at r < i?c = 3.25cr with ip{r) = 
otherwise, where 



$ = 4u 
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(17) 



The MD test runs have been performed at a reduced 
density of n* — 0.845 and a reduced temperature of 
T* = k^T/u = 1.7. For the purpose of comparison, 
the equations of motion were integrated applying also a 



well-established velocity Verlet (VV) algorithm [|T^,^ of 
the second order and its forth-order (VV4) counterpart 
]l8t . Our algorithm we will call now as CPFD (conser- 
vative pure fluid dynamics). A typical maximal value 
for the reduced time step in simulating such a system is 
T-max = T(u/mcr2)i/2 ^ 0.005 ||. All the runs started 
from a well equilibrated configuration and covered an 
identical time interval of t* = t{u/ma'^y^'^ = 50 (cor- 
responding to 10 000 time steps at r* = 0.005). 

It is worth mentioning that the explicit VV integrator 
propagates the phase variables according to the relations 

r,it + t)= r,(<) + v,(t)T + Ut)— + C)(r3), 

Zm 

v,(t + r) = v,(t) + [i,{t) + f,(t + r)] + 0{t^). 

This propagation can be presented as 

{r, (t + r) , V, (i + r) } = D (i , r) { r, (t) , v, {t) }+0{t^), 

where D(t,T) denotes the evolutionary operator. The 
VV4 algorithm deals (similarly to the ED4 scheme) with 
the five stages propagation 

5 

{r,(t-fT),v,(t+T)} = []D(t,aT){r,(0,v,(t)}+O(T5), 

k=l 

where the coefficients are: = ^2 = = ^5 = ^ = 
1/(4 - 4I/3), and 6 = 1- 46 The VV approach needs 
only in one force evaluation (the most time-consuming 
part of the calculations) per time step, pvv = 1, while 
Pvv4 = 5. The CPFD algorithm requires two force eval- 
uation per iteration within the time step, Pcpfd — 2. 

The averaged total energy fluctuations £ = {{E{t) — 
E{0))'^y/'^/{uN) obtained within the CPFD integration 
at the time step t*^^^ = 0.005 and the numbers of it- 
erations of I = 2,3,4, and 8 are plotted in Fig. 3 as a 
function of the reduced processor time Ip (where in this 
case p = pcpfd) needed to perform the run of the men- 
tioned above length t* = 50. The fluctuations identified 
during the integration at the time steps t* = us- 
ing the VV algorithm with Z = 1, 2, 4, 8 and 16 as well as 
the VV4 algorithm with / = 0.5, 1, 2, 4 are also included 
in this figure. The processor time spent to carry out the 
VV run of the length t* = 50 at r* = 0.005 is assumed 
to be equal to unity in our dimensionless presentation 
Ip (where p = 1,2 and 5 for the VV, CPFD, and VV4 
integrators). 

The LJ energy fluctuations damp with increasing I like 
- - r^, and - cxp(-2.40 within the VV, VV4, and 
CPFD integrations, respectively. Up to three intersection 
points corresponding to the VV-VV4, VV4-CPFD, and 
VV4-CPFD curves with the energy conservation levels 
of £1 ~ 6 ■ 10-^ £2 10-^ and £"3 ~ 3 • 10"^ can be 
observed in Fig. 3. So that the usual VV algorithm is 
recommended to be used when the precision £ of energy 
conservation plays not so important role in the compu- 
tations, namely, when £ > £i. The calculation with the 
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help of the VV4 integrator appears to be most computa- 
tionally efficient in the intermediate regime £^ < £ < £i. 
Finally, when a very accurate conservation, £ < £3, is re- 
quired, the best choice is to apply the CPFD algorithm 
because then it becomes to be most economical. 
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FIG. 3. The averaged total energy fluctuations £ as a func- 
tion of the reduced processor time Ip needed for the simula- 
tions of a Lennard-Jones liquid within the CPFD (bold solid 
curve), usual velocity Verlet (solid curve, VV), and fourth- 
order velocity Verlet (dashed curve, VV4) algorithms. 

The CPFD approach can also be used for the predic- 
tion of dynamical phenomena in other many-body collec- 
tions (such as the solar system, for instance) and treated 
as an efficient numerical solver of first-order differential 
equations. The most notorious example (which allows to 
be analyzed analytically) is the equation (Px/dt'^ = x = 
—X describing dynamics of a simple harmonic oscillator. 
This equation reduces to a system of two first-order dif- 
ferential equations, x = v and v = —x, which in turn can 
be reproduced from the first two lines of general equa- 
tion (2) putting formally Vi = x, Vi = v, (p = x'^/2, 
rrii = 1 and J = 0. Then in view of Eq. (14), the time 
propagation reads x{t + t) = x{t) + T[v{t) + v{t + t)]/2 
and v{t + t) = v{t) - T[x{t) + x{t + t)]. The last two 
relations can be solved explicitly, and the result for the 
conservative numerical trajectories is 



x{t + t) = 



;(t)(l-rV4)+t;(t)r 
l + r2/4 

x{t)+v{t)T/2 



V{t + T)=v{t)-T ^^^^^^ , 

whereas the VV solutions are 

x{t + t)= x{t){l - rV2) + v{t)T, 

v{t + t)= v{t) - T[x{t){l - T^/A) + v{t)T/2] . 



Choosing the initial conditions .x(0) = and i(0) = 
f (0) = 1, the above two types of numerical trajectories 
can be compared between themselves and with respect to 
the exact solution x{t) = sin(t) and x(t) = v(t) = cos{t). 
The result of comparison for x{t) is presented in Fig. 4 
at a typical time step of t = 0.05T, where T = 2tt de- 
notes the period of the oscillations. As can be seen easily, 
the conservative solution leads to a better reproduction 
of the original dependence than the VV trajectory, de- 
spite the both CPFD and VV approaches are valid to 
the same 0{t^) order in truncation errors. Therefore, 
additional cancellations of the truncation uncertainties 
are possible due to the exact preservation of the integral 
of motion 2E = x^ + x"^ = 1 along the CPFD trajectory 
(note that maximal VV deviations in E consist about 
20% at T = 0.05r). Similar cancellations of the trun- 
cation uncertainties in microscopic solutions within our 
conservative approach should be expected for other sys- 
tems of differential equations, in particular, for spin and 
pure liquid dynamics. 




FIG. 4. Numerical solutions to the differential equation 
x{t) = —X obtained during our new (solid curve) and velocity- 
Verlet (dashed curve) integrations at the time step t = 0.05T 
with T = 47r. The exact result x{t) = sin(a;) is shown as open 
circles. 



VI. CONCLUDING REMARKS 



One of the most fundamental characteristics in physics 
are the conservation laws. Therefore, it is desirable that 
the numerical methods in computational physics obey 
these laws. Unfortunately, the most popular algorithms, 
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such as predictor-corrector, Runge-Kutta, Verlet, decom- 
position Suzuki- Trotter, etc., being applied to the non- 
Unear many-body problem, do not preserve fundamental 
physical invariants, like energy and angular momentum, 
when these are inherent in the description. 

In the present paper we have tried to remedy such a 
situation and formulated a novel completely conserva- 
tive approach for numerical integration of the equations 
of motion in classical systems. The approach is general 
enough to be used for a wide class of systems such as spin 
and pure liquids, collections of charged particles, etc. It 
can also be considered for the prediction of other phe- 
nomena in physics, astrophysics, chemistry and biology, 
whenever the numerical solutions to systems of differen- 
tial equations are necessary. 

Our main attention in this study was concentrated on 
dynamics of spin liquid models in which additional ef- 
fects with respect to pure liquids are possible because of 
the energy exchange between spin and liquid subsystems 
ip^-pTf. As a result, a new second-order MD algorithm 
(called as CSFD) has been consequently derived within 
the above presented approach. Its greatest advantage is 
that all the integrals of motion existing in the system, 
namely, the total energy, linear and angular momenta, 
individual spin lengths and total magnetization are con- 
served independently on the size of the time step. It is 
worth emphasizing that such a complete conservation has 
been achieved intrinsically, i.e., without the introduction 
of any artificial external forces or numerical constraints. 
Moreover, the resulting algorithm maintains the time re- 
versibility property inherent in the basic equations. This 
is also important for long-duration MD observations be- 
cause the stability of an algorithm is closely connected 
with its time reversibility |23| ]. 

The presented algorithm is implicit, i.e., it requires 
iterative solutions. Thus, when a high precision in con- 
servation is not needed, the CSFD scheme may be less 
efficient in practice than explicit decomposition methods 
1^,0. We have shown, however, on the basis of ac- 
tual simulation of a Heisenberg fluid model that when 
the total energy and magnetization must be reproduced 
precisely, the CSFD algorithm may be in order or even 
more faster than the decomposition integrators. Another 
important feature of the conservative method is that ad- 
ditional cancellations of the truncation uncertainties are 
possible in microscopic solutions due to the exact preser- 
vation of the macroscopically observable integrals of mo- 
tion, as was demonstrated analytically on a simple ex- 
ample of the harmonic oscillator. 

For a particular case when the spin subsystem is ab- 
sent, the CSFD algorithm reduces to a so-called CPFD 
integrator. While this work has been done we have 
learned that this integrator is equivalent, in fact, to that 
developed independently by Greenspan |^ as well as 
Gonzalez and Simo |2^]. These authors, however, con- 
sidered the integration in respect of applying it to me- 
chanical systems when the number of particles is not very 
large. Here we have shown within the LJ model that the 



CPFD integrator can be used with equal success in MD 
simulations of pure liquids. 

The approach presented can be adapted to many- 
component systems, optimized to a multiple time step- 
ping integration and extended to higher-order versions. 
These and other related problems will be the subject of 
a separate investigation. 
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